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A calculation technique in the context of the self-energy functional approach (SFA) and its local 
form, the dynamical impurity approach (DIA)-, will be proposed. This technique allows for a precise 
calculation of the derivatives of the grand potential functional used in the search for a stationary 
point. To make a closer comparison of the DIA with the dynamical mean-field theory (DMFT)-, 
and to demonstrate the proposed technique, we calculated paramagnetic U-T phase diagram of the 
Hubbard model at half-filling, which exhibits metal-insulator transition. 

PACS numbers: 71.15.-m, 71.30.-fh 
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I. INTRODUCTION 

Variational determination of the grand potential func- 
tional has a long history. The seminal paper, made in 
the 1960s^, provides us with the expression for the grand 
potential functional which reads 

n[G] = $[G] + Trln(-G) - Tt{{Go^ - G-^)G), (1) 

where Go is the free Green's function, G is the full 
Green's function and $[G] is the Luttinger-Ward hmc- 
tional expressed by a sum over all connected skeleton 
diagrams. fl[G] is stationary for the physical Green's 
function G. The main problem which arises in that ap- 
proach is a formidable skeleton perturbation expansion. 
Approaches which used Eq. J^l as a starting point would 
approximate the Luttinger-Ward functional by an incom- 
plete diagram series or by composing ^'[G] from only a 
few lowest order diagrams^. Recently it has been pro- 
posed to regard the grand potential functional as a func- 
tional of the self-energy with the Legendre transform of 
the Luttinger-Ward functional calculated from a simpli- 
fied Hamiltonian (which we will refer to as a reference 
system) and the self-energy of the original system ap- 
proximated by the one of the reference systemi. Mathe- 
matically expressed, the proposal states that 



r!t[S(t')] = r!t'[S(t') 



Trln(-(Go-i 



Trln(-(G'-' 



-s(t'))-i) 

S(t'))-'), (2) 



where S[t'] is the exact self-energy of the reference sys- 
tem, Gg is the free Green's function of the same and Gq 
is the free Green's function of the original system. The 
reference Hamiltonian is chosen such that it is simpler 
than the original one, while a systematic increase in its 
size would reproduce the original system. In the case of 
the Hubbard mode&Sti , correlated sites in the reference 
system correspond to those from the original one. We 
can add uncorrelated sites to correlated sites to mimic 
the rest of the latticoi*^. The grand potential functional 
is stationary with respect to the variation of the exact 
self-energy. The same stationarity condition will be im- 
posed to the approximate one defined by Eq. Q with 



the reference system being simpler than the original one. 
The simplest form of the reference system is obtained 
if we restrict our attention to the local self-energy, in 
which case we speak of the dynamical impurity approach 
(DIA). In general, including reference systems with non- 
local self-energies, the method is called self-energy func- 
tional approach. In this account we will concentrate on 
the calculation performed in the context of DIA, even 
though, only slight changes are necessary for a more gen- 
eral, SFA approach. 

In the seminal paperi it was demonstrated on the ex- 
ample of the Mott-Hubbard transition how even a rudi- 
mentary reference system can give a good account of 
the physics obtained by the numerically more expensive 
DMFT approach. We are going to address the question 
of how the extension of the reference system improves the 
results. Since DIA is equivalent to DMFT only when the 
number of uncorrelated sites goes to infinity, it will be 
interesting to compare the convergence of the results ob- 
tained in DIA to those obtained in DMFT as the number 
of uncorrelated sites in the reference system is increased. 
Contrary to the mentioned rudimentary case with only 
one variational variable, the set of variational parameters 
in the reference system with six atoms has five compo- 
nents for the half-filled case. To be able to determine bor- 
ders of the phase diagram precisely, we calculate deriva- 
tives of the grand potential functional analytically. This 
procedure will be sketched in Section |n] together with 
additional calculations given in the appendices. In the 
Section lTTTl we present the results for the Hubbard model. 



II. DETERMINATION OF THE STATIONARY 
POINT 

The main numerical burden of the SFA/DIA method is 
the determination of the self-energy. That part consists 
of the diagonalization of the reference system, calcula- 
tion of the Green's function from the Lehmann repre- 
sentation and, using Dyson's equation, determination of 
the self-energy. The only numerically problematic part 
is the one of the calculation of the self-energy from the 
Green's function (see Appendix IB|| . With the obtained 
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self-energy we can evaluate trace terms along the lines al- 
ready discussed in some earlier accounts^. That approach 
would demand numerical calculation of derivatives for the 
determination of the stationary points of the grand po- 
tential functional. The numerical load for a reliable and 
precise calculation of derivatives can be larga^. We will 
show how the additional information from the diagonal- 
izcd reference system can be used for their determination. 
The biggest advantage of the method can be seen in the 
coexistence region of the phase diagram. 

Before embarking on the calculation of derivatives of 
the trace terms, we are going to sketch a few steps in the 
calculation of the trace terms already discussed^. Using 
a general form of the Green's function G, the trace term 
can be written as 



Trln(-G) = ^— ^trln 



-1 



M - t - S(ia;„) ■ 



(3) 



where tr denotes trace over all quantum states. A 
function of a Hermitian operator is defined as f{A) = 
Uf{a)U^ where a is the diagonalized A matrix and U is 
the unitary transformation which performs it. If r]k{iujn) 
are eigenvalues of the operator t -I- S(itt'„), then 



In 



(4) 



Summation over k encompasses all eigensolutions of the 
operator t + H{iuj). In order to get over from the summa- 
tion over Matsubara frequencies to the integration in the 
complex plane, a standard trick is applied. It consists 
in the introduction of a function in the complex plane 
which has poles at z = iujn and the subsequent integra- 
tion around the contour which includes the imaginary 
axis 



Trln(-G) 



■In- 



-1 



^ Jc 27ri ef^^ + 1 z + - rik(z) 



(5) 



The integration contour in Eq. ^ can be mapped to the 
integration contour enclosing the real axis by the lines 
infinitesimally above and below it {u)±i5). Using analytic 
properties of the integrand, i.e., 77fe(w — i5) = ri^{uj + iS), 
and retaining only the largest contributions in 0"*" and S 
(both are infinitesimally small positive numbers), we can 
express Trln(— G) as 



E 



-duj 



/(ti;)Imln 



uj + iO+ + fi- 77fe(cj + iO+) 



(6) 



where /(w) = l/(e^'^ + 1) is the Fermi-Dirac function. 
For a function F{w + iO^) that has the property ImF(a;-|- 
iO'^) < (that is the case for the diagonal Green's func- 
tion), there follows the equality Imln(— F(a; -I- 10+)) = 
ttQ{F{u;)) = 7re(l/F(w)), with 6 being the step func- 
tion. Due to the finite number of poles in the reference 
system, singular points don't give a finite contribution to 
the integral. The resulting expression for the trace term 
reads 



Trln(-G) = - 



dujf{uj)Q {uj + ~ r]k{io)) . (7) 



This expression is our starting point for calculating 
derivatives with respect to the parameters of the refer- 
ence system. As in the calculation of the grand potential 
functional^, we want to obtain the integrand as a sum of 
the Dirac 5-functions, which is then trivially integrated. 
For this purpose we just have to exchange integration 
over the real axis and derivative with respect to the ref- 
erence system parameter A. Due to the discontinuity of 
the argument of the step function in Eq. |(7J), we obtain 
two sums of the (5-functions, one evaluated in the zeros of 
uj + ^ — ryfc(S(t'), u) and the other evaluated in the poles 

0f7?fc(S(t'),C^). 



d_ 
dX 



Trln(-G[E(t')]) = -^^ j dc./(u;)e (^ + /i - 7?fc(S(t'), c.)) 

/oc 
-oo 



5[uj + i-i- r]k[uj)) 5 



cj + /i-?7A;(S(t'),cj; 



k.i 



= J2 / dijfiu 



97?fc(S(t'),c.) S{co-J,^^) 



dX 



duj 



d 


(w+Ai-r,fc{S{t').t^)) 






dX 







L/K.) dx 



duj 



(8) 



dX 


1 - 


a»;fc(S(t'),c^) 






duj 


UJ—LO ' 



(9) 



,(10) 
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The second term in Eq. (|10|l should be thought of as an 

evaluation of the derivatives for lo — > since each of 

(v) 

the derivatives diverges for uj = Uj ^. . Using the notation 



of Eqs. (|B13I IbT4|) . but with coefRcients being matrices 
instead of scalars, we can obtain the final expression for 
the derivative of the trace term. It reads 



— Trln(-G[E(t')]) 



E 



9%(S(t') 



8r)fc(S(t'),i^) 



(11) 



(1) 5j 

where \'^k.j) is an eigenstate of the matrix . Equa- 

tion is completely general, i.e., it can also be applied 
to the case of non-local self-energy. From now on, we are 
going to apply it only to the case of the local self-energy. 

For obtaining the grand potential functional we have 
to evaluate two trace terms: one for the Green's function 
of the original system and one for the Green's function 
of the reference system [see Eq. (^]. Reference systems 
with the local self-energy can have any form (chainlike, 
starlike etc.), because reference systems with different 
configuration of bath sites can be mapped to each other 
so that the physical Green's function doesn't change. We 
have chosen the star-like reference system for its sim- 
plicity, with the Hamiltonian given as (the paramagnetic 
state is assumed) 

+ ^(ea,z - ^J)a\aala + ^Vi{c\^^ai^„ + a\^^ci^^). (12) 



The trace of the full Green's function of the reference 
system can be decomposed in a number of trace terms 
with the scalar Green's functions 

tr' ln(-G'(iw)) 

No. 

= J2 H-G[A^u;)) + H-GlA^u;)), (13) 



cr CT k=2 

with G'i^{iLUn) = ^/{i^n + M — Efc) and 

1 



A* - ei - E 



k—2 iujn-\-^—^k 



T.{iuj„) 



-(14) 



Using the notation of Eq. (|10|1 . the following identifica- 
tions can be done: 



r]i{iujn) = ei + E 



k=2 



+ S](jc^„) (15) 
Vk>i = Efc (16) 



The original Green's function is diagonalized by the 
Fourier transformation under assumption of the local 



self-energy and translation invariance of the system. In 
this case rjk,a{i^) = e(^) + S(w), where e(fc) is the Fourier 
transformed one-particle term and k denotes a wave num- 
ber. 

Stationary point in the parameter space t' assures us 
that we have found the self-energy which is physical. 
Since it is possible to find more than one solution, ther- 
modynamically stable are only those characterized by the 
minimum of the grand potential^. 



III. METAL-INSULATOR TRANSITION IN THE 
HUBBARD MODEL 

As a demonstration of the above numerical proce- 
dure, we determined the U-T phase diagram of the Hub- 
bard model in the paramagnetic state. In addition to 
the results for two-site DIA which have shown qualita- 
tive agreement with the full DMFT calculations, we will 
demonstrate that richer reference systems allow us to 
make quantitative comparisons with the existing DMFT 
calculationsiSiiiiiS. 

In contrast to the DMFT, where the local Green's func- 
tion of the reference system is identified with the one of 
the original system, we have to distinguish the two in 
the DIA formalism. The difference between DMFT and 
DIA procedures is best seen from the stationarity con- 
dition of the grand potential functional which results in 
the identity 



EE 



1 



Go 1 - S(f ) 



G' 



0.(17) 



5n.m^n)- This 



The self-energy in DIA is local (I]„^„ 
implies that the stationarity is automatically fulfilled if 
the Green's function of the impurity site in the reference 
system is identical to the local Green's function of the 
original system. The identity of the two Green's hmc- 
tions in DIA formalism is expected only when the num- 
ber of additional bath sites is large. In general, the local 
Green's function of the original system and the impurity 
Green's function of the reference system are related to 
each other by the self-energy which they share, but oth- 
erwise they are different. Thus to satisfy the stationarity 
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condition we also have to take into account the deriva- 
tive of the self-energy with respect to a parameter of the 
reference system appearing in Eq. H17() . In order to fa- 
cilitate the search for a stationary point, wc reduced the 
parameter space using particle-hole transformation which 
is a symmetry transformation of the Hubbard model at 
half-filling. This transformation is defined by the idcn- 

TABLE I; Constraints on the parameters of the reference 
system imposed by the particle-hole symmetry at half-filling 
(ei = Ec, ^i>2 = ^a,i)- N is the size of the reference system. 

N Variables and their relations 
'2 ^yV2 

3 a) V2,V3, £2 = £3 = 

b) Vi = V2,e3 =2/1-62 

4 a) V2,V3,V4,e2 = £3 = £4 = ^ 

b) V2,V4 = 1^3, £4 = 2^ — £3,62 = M 

5 a) V2, V3, 14, V5,e2 = £3 = (;4 = £5 = M 

b) V2, Vs.Vs = V4,e2 = £3 = iJ^,^5 =2/1-64 

c) V3 = V2,V5 = Vi.ez = 2fi- £2,65 = 2/i - £4 

6 a) V2, V3, V4, V5, V6,£2 = £3 = £4 = £5 = £6 = M 

b) V2, Va, V4, Ve = V5,£2 = £3 = £4 = M.f^e = 2/i - £5 

c) V2, V4 = V3, Ve = V5,£2 = /i,£4 = 2/1 - £3, £6 = 2/i — £5 



tity Pci^aP^^ = e'-'^^cl^, where 0i's are chosen such that 
4>i — 4>j = TT and i,j are site indices. We can choose 
(j) — on the sites with even Manhattan distance and 
(j) = TT for the remaining sites. The same transforma- 
tion is performed on the impurity site. For the bath sites 
we have chosen such a particle-hole transformation that 
the resulting number of variational parameters is maxi- 
mal. For the case of the impurity site with cj) = 0, bath 
sites should be transformed according to the expression 
Pai^^p-^ = -a|^. The transformed if^/AA/ Hamilto- 
nian is then 

PHP-^ =^(/i-e, -[/)4c, + C/ntni 
+2 ^(e/ - m) + C/ + 2(e, - /i) + ^(/i - ea,/)^.^,,. 

Under the condition that the transformed Hamiltonian 
and the original Hamiltonian of the single impurity An- 
derson model (SIAM) are the samei^, we can put con- 
straints on the parameters of the reference system. Since 
/i = J7/2, we find immediately that tc = 0. Bath states 
transform according to the equations: 

M - Ca,; = Ea,;' - Vi = Vv. (19) 

Combinations following from Eq. H19|) are listed in Ta- 
ble |l] for different sizes of the reference system. It can be 
shown that the effective size of the reference system is 
determined by the number of different e's, on-site ener- 
gies of the bath orbitals. Namely, for the orbitals with e's 
equal to the chemical potential it is possible to perform 



a unitary transformation that transforms uncorrclated 
sites with ei = fj, into one site coupled to the impurity 
site and the rest of the uncorrelated sites uncoupled from 
the impurity site. For example, if we had two uncorre- 
lated sites in the reference system with £2/3 = M then it 
can be shown that the Hamiltonian can be mapped to 
the one with the hopping amplitudes V2 = a/Vj^ -I- V^, 
1/3 = 0. Thus the effective reference system has been 
reduced to a system with only one site. In the case of 
more bath sites with a = fx we can iterate the transfor- 
mation for two sites and show that the resulting effective 

amplitude is V2 — \l J2iL2 j ^>2 ~ 0. The solutions in 
this work correspond to the cases 2a, 4b and 6c from Ta- 
blenievaluated for the semicircular free-particle density of 
states. In Fig. ^ we show the calculated quasi-particle 




FIG. 1: Quasi-particle weight at T=0 calculated in DIA for 
different sizes of the reference system (N=2,4,6). Curves 
for (N=2,4) have already been presented-. For compari- 
son, DMFT results obtained using numerical renormalization 
group(NRG)-- and exact diagonalization(ED)- as impurity 
solvers are included. 

weight Z = 1/(1 - aReE(w + iO+)/dLu)\^=o. Regarding 
the interval of U values as a whole we can notice that the 
convergence of Z with increasing system size is uniform. 
That is not true in the vicinity of the phase transition. 
Even though the convergence is not uniform in that re- 
gion, it seems to be fast. The comparison with DMFT 
calculation using exact diagonalization as impurity solver 
shows a close connection of the two methods. The phase 
diagram of the paramagnetic state of the Hubbard model 
in the U-T plane is shown in Fig. 0. Three different re- 
gions in temperatures below the critical point, denoted by 
an empty circle, can be distinguished: (a) metallic phase 
for small values of U, (b) insulating phase for large U, (c) 
coexistence of both phases in a triangle-like shape for the 
intermediate Coulomb interaction. Whereas the transi- 
tion below the critical point is of the first order, we find 
the crossover from the metallic to the insulating solution 
in temperatures above the critical point. Convergence of 
the phase diagram for different sizes of the reference sys- 
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tern is shown in the inset of Fig. On the metalhc side 




FIG. 2: U-T phase diagram of the paramagnetic state of the 
Hubbard model. Ud curve denotes the left border of the re- 
gion with the insulating solution. Uc2 is the right border of 
the region with the metallic solution. Coexistence region has 
both solutions. The insulating solution is stable to the right 
of the curve Uc, whereas the metallic solution is stable on its 
left. The empty circle indicates the critical point separating 
low and high temperature regions with the first and the sec- 
ond order phase transitions between metallic and insulating 
solution respectively. Uci and Uc2 curves for different sizes of 
the reference system(N=2,4,6) are shown in the inset. 

of the metal-insulator transition in the DMFT formahsm, 
the central role is played by a three peak structure in the 
spectral function, the middle peak corresponding to the 
Abrikosov-Suhl resonance in the impurity model and two 
Hubbard bands. That distribution of the spectral weight 
together with Tabled can explain the convergence trends 
in the solutions found for different sizes of the reference 
system. For N = 2 we have only one pole in the inverse 
of the free Green's function of the reference system at 
the Fermi level. From the DMFT equation for the semi- 
circular density of states GQ"^(aj„) — iujn + A* ~ t^G„{u}n) 
it would follow that the local Green's function has only 
one pole on the Fermi level if the equation holds for each 
number of bath sites. In the DIA formalism the connec- 
tion between the on-site Green's function and the inverse 
of the free Green's function in the reference system of a 
finite size is more involved, but we believe that rapid con- 
vergence of the DIA results to those obtained in DMFT 
allows us to use the DMFT equation for the argumen- 
tation. This means that, with N=:2, wc cannot properly 
account for the Hubbard bands central to the insulat- 
ing phase. The Uci curve is thus substantially under- 
estimating the extension of the insulating phase region 
in comparison with N = {4, 6} results. The same rea- 
soning explains why N = {3,5} reference systems show 
no improvement of the metallic solution with respect to 
the reference systems with a smaller but even number of 
sites. A reference system with an odd number of sites. 



due to the particle-hole symmetry, has no pole in the in- 
verse of the free Green's function at the Fermi level or it 
has two poles at the Fermi level which can be mapped to 
one. For = 4 we can account for the side bands, and 
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FIG. 3: Shifted grand potential functional in the direction of 
the parameters of the reference system around the insulating 
solution for T=0.0I6 and U=4.77. Grand potential functional 
is given by n = fis^ -1-2.50717. Schematic configuration of the 
reference system with the corresponding parameters is shown 
in the upper left corner. 

A = 6 brings only a slight change to the phase diagram. 
Boundaries in the phase diagram are defined by the dis- 
appearance of a stationary point in the parameter space 
for either metallic or insulating solution. An example 
of the parameter space is shown in Fig. The calcu- 
lation was done for six atoms. At half- filling and in the 
paramagnetic state there are five independent parameters 
(upper left corner of the figure). The parameter space is 
shown in the region around the stationary point for the 
insulating phase and U=4.77, T=0.016. Insulating so- 
lution disappears when the stationarity condition is not 
any more fulfilled in the direction V2. As already argued, 
V2 is related to the weight of the Green's function of the 
original system at the Fermi level. Another comparison 
of the calculations done in the DMFT framework with 
the DIA calculation is shown in Fig. Q for the whole 
phase diagram. As already noticed for T=0, the DIA 
calculation shows strong resemblance to the DMFT-ED 
results. 

In conclusion, wc showed how additional information 
from the reference system can be used to increase the 
precision of the calculation in the context of DIA. Com- 
parisons with the results obtained in DMFT demonstrate 
close connection between the DIA and DMFT-ED pro- 
cedures. The advantage of the DIA formulation over the 
DMFT-ED formulation is particularly clear for the case 
when the number of bath sites is not sufficiently large 
to reproduce the local Green's function of the original 
systemi. Already for A = 6, DIA and DMFT-ED give 
almost the same result for the metal-insulator transition 
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— DIA-N=6 
□ oDMFT-ED 
* ODMFT-NRG 
A ADMFT-QMC 




FIG. 4: Metal-insulator transition in the half-filled Hub- 
bard model. The solid line shows the result from DIA 
with a reference system of the size N=6. DMFT results 
with ED(squares)--,NRG(diamonds)— and quantum Monte 
Carlo(triangles)— as impurity solver are shown for compari- 



calculated by applying the derivative to its Lehmann rep- 
resentation. For each pair of indices {i,n) there are two 
terms contributing to the Green's function, where the 
first term is the projection of the state created out of 
the N-particle ground state upon addition of a particle, 
on the eigenstates of the (N-|~l)-particle system. The 
second term is the projection of the ground state with 
added hole. We can get the second term by performing 
the following replacements in the first term: c^, 
cl ^ Cb, and by exchanging energies {Eq <^ En) in the 
denominator. The derivative of the first term is 



d_ 
dX 



(n|ct|nW)(„|ct|n«) 



LJ-iEn- Eo) 

I ti 
-(n\cl\n^ 



g I 



{n\cl\ng 



dX 



Lo ~ (S„ ™ Ea) 
(n|4|4'))(n|c^|4'^) fdEn dEo 
[lj - {En - EoW \ dX d\ 



(A2) 



in the paramagnetic state of the Hubbard model at half- 
filling. 

I would like to thank Prof. C. Gros for his help in a part 
of this work, M. Potthoff for kindly providing additional 
information on the SFA prior to the publication, and W. 
Krauth for a stimulating discussion. 



In the term with lo — [En — Eq) m the denominator, the 
second part in the numerator is obtained from the first 
The remaining derivatives are: 



by changing to cj. 



d{n\cl\4^) 



dX 



dEn 
dX 



-^ciy;>) + {n\cl^^, 

dEp 
dX 



(A3) 



APPENDIX A: 



DERIVATIVE OF THE GREEN'S 
FUNCTION 



In Eq. O and all subsequent ones we evaluate the real 
part of the self-energy and the Green's function on the 
real axis (poles are excluded) . Thus we will be interested 
only in the derivatives of the real parts of those func- 
tions. Since Green's functions are calculated from their 
Lehmann representation, 



"l4l4'^)('^|Cbl»g' 
UJ-{En- Eo) 



+ 



(n|cb|4*^)(n|ca|4'^) 



OJ 



{Eo — En) 



(Al) 



where iV,; is the degeneracy of the ground state, 
are degenerate ground-state eigenfunctions and Eq is the 
ground-state energy, the derivative with respect to lo is 
straightforward to evaluate. Summation over n extends 
over all eigenstates. In Eq. (|A1|I we used the fact that 
the matrix representation of the Hamiltonian and the 
corresponding eigenfunctions will be real, and organized 
all matrix elements so that the ground states are in the 
ket state. 

The derivative of the Green's function with respect to 
a parameter of the reference system A G {ei,T^} can be 



In the end we have the problem of determining deriva- 
tives of eigenvectors and eigenvalues with respect to the 
parameter A. It is well known that the terms up to the 
second order in the Rayleigh-Schrodinger perturbation 
theory correspond to the order of perturbation A. We 
can thus use the perturbation to determine derivatives of 
eigenvalues and eigenvectors. Special attention should be 
paid to the case when some eigenstates are degenerate. 
If we denote the perturbation as H' = XV, then, in the 
perturbation theorjii^, we have to use states from the 
subspace which diagonalize the operator V. Assuming 
that the above requirement is fulfilled, it follows that 



d{n\ci\4^) 
dX 



{m\V\ng^''){n\cl\m) 
Eo — Em 



E 

m 

^ {n'\ci\4^){n'\V\n) ^^^^ 



En — En' 



and 



dEn 



= {n\V\n) 



^^°=(n«|y|4)) (A5) 



ax ^ ' ' " dx 

After some rearrangements, the derivative of the Green's 
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function has the form 
dRcGab{uj,\) 1 



dX 



1 ab 



EE 



B: 



ab 



uj - {En -Ea) Lu- {Eo - E,,) 



b 



1.71 



[lu - (S„ - EoW - (^0 - ^n)]2 j 
with weights of the first order poles given by 

= E |^((-'l4l-«)("|ci|n«) 

n 

+ (n|ct|nW)(n'|ct|^W)) 

+ E ^^^^(("l4HH4k 

+ (n|ct|n«)(n|ct|m) 



(A6) 



(A7) 



exchange of the operators. The second order poles have 
the weights CfJ'j equal to 



(n|ct|nW)(n|ct|n«) [{n\V\n) - (nW|y|nW)) (A8) 

and DfX = -Cliicl^Ca,cl^Cb). 

Along the same lines we can calculate the derivative 
of the Green's function for finite temperatures. For the 
Green's function ReG'ah(w) being given by the Lehmann 
representation 



E 



1 (^-PE„' + g-''^") ("'l^°l")("l^bl"') 
— {En — En') 



(A9) 



the derivative ox'' given by the function 



E 



Aab 
n.n' 



iDab 



(AlO) 



and Bf^ = A°^j(cjj 4^ Ca,cl Cb), where <^ means the with AJ^''^, equal to 



n|4|n')(n|4|m) 



E 



I 



{l\V\n) 

En — El 



{l\ciW){n\clW) + {n\ci\n'){l\cl\n'))+ ^ 



i\V\n') 



En' - E, 



1 

Z 



P 



dEn 1 dZ 



d\ Z d\ 



a£:„' 1 dZ 



d\ Z dX 



(n|4|m) {n\cl\n') 
n\ci\n'){n\cl\n') (All) 



and 

<„, = ^{e-^^'^+e-P^-'){n\cl\n'){n\cl\n') 

X i{n\V\n) - {n'\V\n')) . (A12) 

APPENDIX B: SELF-ENERGY AND ITS 
DERIVATIVES 

Self-energy 11{t',uj) is defined by the Dyson equation 

^{uj) = Ga\u)-G-\uj), (Bl) 

In order to be able to perform high precision calcula- 
tion, it is necessary to be able to do exact subtraction 
of Gq^{u;) and G~^(w). Direct subtraction of the calcu- 
lated inverses of Go{uj) and G{uj) at some point to is not 
a method of choice if we strive for high precision calcu- 
lation. The reason lies in the finite precision arithmetics 



in each computer. Details of this discussion will be given 
at the end of this section. 

An alternative to the subtraction of the evaluated in- 
verses of Green's functions is to subtract them as func- 
tions of the frequency and only then evaluate them in 
the frequency we are interested in. The difference with 
respect to the subtraction of already evaluated functions 
is that in this way we always subtract numbers of or- 
der one, whereas in the direct subtraction the subtracted 
numbers are sometimes large (in the vicinity of a pole) 
but close to each other. In this section we will focus only 
on the DIA case, i.e., local self-energy. Generalization to 
the non-local self-energies is straightforward. 

If we mark the impurity site with index 1, then the 
self-energy on that site has a form 

ReEn(t^) = [Go,ii(^)]"' - [Gn{uj)]-' (B2) 

where, in the case that each bath atom is connected only 
with the impurity site (starlike bath), the inverse of the 
free Green's function is obtained from the equations of 



8 



motion and reads 



Nt 



[ReGo,ii(w)] ^=w + At-ei-^ 



1=2 



Lu + n - ei' 



(B3) 



where e/ is on-site energy, Vi is the hopping amphtude be- 
tween bath sites and the impurity site, Nb is the number 
of atoms in the reference system and /i is the chemical po- 
tential. From the Lehmann representation we know that 
the real part of the thermal one-particle Green's function 
on the real axis can be written in the form 



ReGiiM =^ 



Ng ^(1),G 



(B4) 



where cf^'*^ > and ^, ^1'''" = 1- 

is a number 

of poles and e^'s arc poles of Gu{u!). The inverse of 
RcGii(a') is then 



(1),G 



The coefficients C. 

i(l),dG-i 



{l),dG 

d 



and C, 



(2),dG" 



follow from 



-{uj-ujjf (9RcGii(w) 



[ReGii(a;)]2 d\ 



(1),G- 



-2c: 



aReGii(w) 



p(0) 



G- 



c: 



W-G- 



c. 



dX 



c 



m-G- 



E 

i=l 



i\2 a^RcGii 



lim 



\ dX 



dujdX 

9ReGii(w) 



(B9) 



dX 



(BIO) 



Unfortunately, the above equation for G^^-'^'^*^ turned 
out to be numerically unstable, i.e., in practice we get 
only a few relevant digits correctly. There is another 
form in which we can write it and which is numerically 
stable. Taking the derivative of Eq. ljB5|l with respect 

to A, we see that G, 



(l),dG~ 



Eq. ljB6|l we can show that C 



^(^ji),G ^Q-^ Now using 

(l),dG 



equals to 



[ReGnH]-'=G(°)'^~'+c.+ J] 



Ng-1 ^(i),g- 



(B5) 



where Wj's are zeros of G (poles of G ^). Coefficients 



G 



(l),G-i 



are obtained if we equate the expression in 



Eq. I)B5|I with l/ReGii((jj), where ReGii(a;) is given by 
Eq. l|B4p , and then calculate the limit lo Uj. We get 



Cf^'°~' = lim 



^Ng C,'^''^ 



(B6) 



To have a value associated with cf^''''^ , we once again 



equate the expression in Eq. I|B5|) with l/RcGii(Li;) from 
Eq. (|B4|I in any non-singular point uj, but this time we use 



just calculated G, 



(l),G-i 



coefficients. The procedure for 



the calculation of the derivative of ReG^^ with respect to 
some parameter of the reference system A resembles the 
one used for the calculation of ReG~^. From Eq. ljB5|l we 
see that i9[ReGii(w)]~^/9A can be written in the form 



9[ReGii(a;)]-i 
dX 



dRcGniuj) 



[ReGii(w)]2 

E ^ 

i=l 



dX 

Ng-I ^(2),dG'^ 



(jj — (jJi 



z— 1 ^ ' 



(B8) 



g: 



(1).G- 



2 / a2ReG 



duodX 

From the definition of w^. 



.(Bll) 



dX duj^ 
i.e., zero of the Green's func- 



tion, its derivative can be deduced and equals 



dujj_ 
dX 



= -G, 



(i),G-i cMeG 



dX 



(B12) 



Coefficient G(°)'''^ t hen fo llows from Eq. jHHl) with the 

coefficients Eq. (jBllI IBlOp and Eq. (|B7p if we equate 
them in an arbitrary non-singular point oj. Having the 
inverses of the Green's functions and their derivative ex- 
pressed as functions of frequency allows the representa- 
tion of the self-energy with a sum of poles added to a 
constant term. Causality of S demands that each pole 
has a positive residue. As a consequence, the poles of 
[Go,ii(ti-')]~^ in Eq. I|B4|) are compensated by the poles of 
[Gii(ci;)]~^ from Eq. ljB5|) . The final expression for the 
self-energy can thus be written in the form 



Rcl]ii(w) = G("'^^ + 



E ^ 

i=l 



UJ ■ 



(B13) 



where G^"^^^ 



G 



(1),G- 



ei - G(")-G and G\ 
1 

. Analytic form of the deriva- 



tive of the self-energy '^"^'^^^ ^"-^ is determined by 



(B7) Eqs. (lijl.HlijHllj8ll 



C<(o) 



E 



G 



U! — UJi 



dX 

Ng-1 

-E 

i=l 



G 



(2),dE 



(B14) 
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TABLE II: Illustration of the degradation in precision by the 
evaluation of the self-energy for different frequencies using 
evaluated inverses of the Green's functions in Eq. IBll {2"'' 
column) in comparison with the result obtained from Ea. lBlSl 
(3'^'* column). The calculation is done for the single impurity 
Anderson model with one bath site and for the parameters: 
U=4, = U/2, = 0, T = 0.001, ea = fi, V2 = 0.5. 





numeric 


analytic 


-le-07 


2.00278960194782 


2.00000017777778 


-5e-08 


2.01115762927247 


2.00000008888889 


-le-08 


2.27893444799702 


2.00000001777778 


-4e-09 


3.74339942335428 


2.00000000711111 


-2e-09 


8.97438195502764 


2.00000000355556 


7.7e-23 


1.6761776955e-Hl5 


2 


2e-09 


8.97277702135762 


1.99999999644445 


4e-09 


3.74317856255584 


1.99999999288889 


le-08 


2.27895191107018 


1.99999998222222 


5e-08 


2.01115782492434 


1.99999991111111 


le-07 


2.00278900765716 


1.99999982222222 



with 

^(0),dE^_^_^(0),dG-i^ (B15) 

oX 

Q = -2 2^(5„,^,Vi— - Q ' , (B16) 
I 



The derivative of the self-energy with respect to the fre- 
quency UJ is immediately evaluated from Eqs. (|lj3IIj5() 
and reads 

As commented at the beginning of this section, the above 
scheme is numerically more stable than the direct sub- 
traction of the functions evaluated at a certain frequency 
UJ. This is illustrated in Table Hll With a given precision 
of 8-byte floating point number (about 15 relevant dig- 
its) we see that direct evaluation gets the first relevant 
digit wrong already when the frequency is at a distance 
of « 10~^ from a spurious pole. Instead of having a finite 
value at = 0, numerical calculation diverges. The non- 
existing divergence in the subtraction of the evaluated in- 
verses of Green's functions has been removed in the step 
where we subtracted coefficients of G^^ and . As a 
result we have a small weight of the self-energy. If we had 
not nullified all small weights (small in comparison with 
the strongest pole) , we would get very weak divergence in 
that point, but we know from the evaluation of the trace 
terms that weak poles give a small contribution to the 
grand potential. This allows us to nullify weights smaller 
than the precision determined primarily by the integra- 
tion performed in the Trln(— G^) term. Aside from the 
loss of precision in the numerical algorithm, the appear- 
ance of the singularity at a point where the self-energy 
is supposed to be analytic might cause problems to the 
zero-searching algorithm. 
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